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In the main program listing: 

page 72 Change card number 155 to read 

DO 601 N=1 , NUMNP 



155 



page 73 Remove card preceding card number 180 that reads 
189 CALL SYMSOL (1) 

page 74 Change statement number on card number 203 to read 
1020 FORMAT (6I5,D10.3) 



In SUBROUTINE FORM: 

page 76 Insert a missing card between card numbers 285 
and 287 to read 



CV «• SPHT(MTYPE) * CFUNC(MTYPE,TMEAN) 



286 



In the MAIN Program listing: 

page 72 Change Card Number 138 to read 
TEMP=TEMPB(N ) - T0 



138 



Change the 


+ signs to - signs 


following 


Card Numbers: 


144 

145 


^•(Note: + Q is OUT ) 


146 




147 




149 




151 
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I* INTRODUCTION 



In this chapter the reasons and advantages of the finite element 
analysis method are given. The purpose of this presentation and the ob- 
jectives of the study are discussed. 

A. REASONS FOR THE FINITE ELEMENT ANALYSIS 

With the ever increasing complexity of problems occuring from engi- 
neering situations in design and analysis, there is a rising need for an 
approximate method of solution for these problems. This technique must 
be flexible enough to encompass a large variety of problems and still give 
accurate answers to the simplest of situations. Since most approximate 
methods of analysis of complex engineering problems involve many tedious 
calculations, it would enhance the usefulness of the technique if it were 
capable of being programed for computer solution. This would relieve 
the engineer of the tedious burden of hand calculated solutions and leave 
him more time to devote to the job of engineering. 

The finite element method of analysis is an approximate method of 
analyzing engineering problems which meet these specifications. It has 
a high degree of flexibility and can be readily applied to computer sol- 
utions to obtain rapid solutions to complex problems. 

B 0 OBJECTIVES OF THIS STUDY 

The major objective of this study is a computer program using the 
finite element method of analysis, for the solution of two-dimensional 
unsteady heat conduction problems. In fulfilling this goal there were 
several minor objectives accomplished along the way. The original pro- 
gram was written for the IBM 7094 system. It was converted for operation 
on the IBM/OS 360 computer at the Naval Postgraduate School Computer Center. 
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Involved in this conversion was the changing of control cards and a few 
FORTRAN IV language differences. Also faster random access disks were 
substituted for the magnetic tapes originally used for temporary storage. 
The input and output formats were overhauled to make better use of space 
on the printed output page. More extensive labels were added to the 
printed output in order to ease the use of the program by persons not 
familiar with the finite element technique. 

C. PURPOSE OF THE PRESENTATION 

The purpose of this presentation is to make available a computer 
program using the finite element method of analysis for the solution of 
unsteady heat conduction problems. In doing so an introduction into the 
formulation of the finite element method is given for those not familiar 
with the technique, to give an indication of the power of this type of 
analysis. A description of the components of the program is included to 
enable the user to better understand the logic of the computer solution 
and increase the user's ability to formulate problems to which the pro- 
gram can be applied. Specific instructions on formulation of problems 
and loading of input data is made available in a chapter on user infor- 
mation, Finally, as is customary in works of this type, recommendations 
for future modifications of the program and possible applications are 
included. 
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II. THE FINITE ELEMENT FORMULATION 



In this chapter the finite element method is formulated. The idea 
of a functional of a function and the variation of a functional are intro- 
duced. These principles are developed into an approximate solution of 
the two-dimensional unsteady heat conduction equation using a matrix 
iterative solution. 

A. FORMULATION OF THE VARIATIONAL METHOD 

If a functional of a temperature field can be found such that its 
extremum yields an equation describing unsteady heat conduction. 



and an equation describing the heat flux across the surface of the region, 



then the variation of the functional will yield a set of first order, 
ordinary differential equations in terms of the nodal temperatures. Using 
the method of tensor notation, a comma denotes differentiation with res- 
pect to the following subscript and repeated subscripts imply summation. 

Let IT be a functional of the temperature field T(x,y,z,t) and the 
first time derivative of the temperature field T(x,y,z,t), be defined by 




(2-1) 




( 2 - 2 ) 




(2-3) 



The variation of 77*(T,T) with respect to T (with T held constant) 
is given by 



<$€ 



(2-4) 
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where £ is a small parameter and A is one of a family of functions. 

A is zero on the boundary of the region and arbitrary elsewhere. An 
extremum of the functional TT(T,T) implies that Sica 9 T) equals zero, i.e. 



<fa-(T-efcA,T) 
c£ £ 



e -o 



o 



(2-5) 



First 



TKT + a,t)= 



V 



’J 



+ f c.(J 4- e X)T - t(T A) )iv 



Jn; V cr + a)ds 



(2-6) 



Using the identity 

1 |cr*«*), s n rj <r*t».J]= [(t*- H.jJ 

if the conductivity tensor is symmetric i.e. kij « kji. Then 

fa (r*tX,t) = A>,- kij X] “ kyCT^feX),.. A 

6e J Y ^ 

f^cXT-^X^dv “ l (2-7) 

P ^ 

The volume integral \(X Rt;A) lT <Jy can be transformed into a surface 

( 2 - 8 ) 

s J 






^ / 9 



integral. 



f ktj X),j dv = 
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When (2-7) is evaluated at £ « 0 

£tt(T + 61/t) -KcjTi-A. 

or 

<$Tr(T-vfrA,t) = J (-kijT,^ +fcT-i')Xd. 
f J Vtaj IjjXdS - J n^ ; X<Js 

The extremum of TT O cT)r educes to 

HcjT, iC = (€t-% 



(2-9) 



( 2 - 10 ) 



in the region V. Also on the surface S 

^sjTvj = nvV 



( 2 - 11 ) 



These two equations (2-10, 2-11) are the unsteady heat conduction equation 
and the boundary flux equation. This verifies that a function T which 
yields an extremum of the functional defined by (2-3) satisfies both the 
unsteady heat conduction equation in the region V and the boundary flux 
equation on the surface S (jQ*. 

B. FORMULATION OF THE FINITE ELEMENT PROBLEM 

# 

If the choice of T and T is restricted so that they are represented 
by certain constants which are obtained in their formulation, the func- 
tional TT becomes a real-valued function. In this case IT shall be 
TfdtrMt}) where £t] and ^t} are vectors of nodal point values 
of The extremum of must now be found which 



^Numbers in brackets refer to similarly numbered items in Bibli- 



ography. 
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implies 



_ O (2-12) 

cST; 

From this point on the finite element formulation will be considered 
in two dimensions only. The region V will be divided into a number plane 
triangular elements. The functions T and T can be uniquely described 
within each of these elements by the values of the function at each of 
the nodal points of the element T^, Tj , T m , and linear function of the 
coordinates of the element. 

Let the temperature field in an element be given by 



TU,v,t) = <0 (x,y)>[a][T ; ®} 



( 2 - 13 ) 



and the time rate of temperature change given by 



TO<,v,t)=<«S(x,v)>[A](t f (t)] 



( 2 - 14 ) 



<0> is a vector which specifies the spacial approximation of T and T. 
The matrix M is a constant and is defined by the above relationships. 
Its exact value and derivation will be discussed in detail in section II, 
E. 

The temperature gradient T,^ ( f| « x,y) can be expressed in terms of 
the nodal temperatures by 






i T ,y 



v ZZ. 






( 2 - 15 ) 



In the following development the subscript ^ , denoting nodal point 
values, will be dropped and will be assumed. 
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Writing the functional IT in terms of nodal point values, 

-Vrf[A]Wi)dv (2-ii 

Taking the variation of IT with respect to ^t| and setting it equal to 
zero 



It], m) = -(£*} =o (2 . l7) 



where 



M= l wVj'WMWJv 



[C] - \ 



V 



[11 = \ iwVdv 



\f} = 



(2-18) 



(2-19) 



( 2 - 20 ) 



( 2 - 21 ) 



C. BOUNDARY CONDITIONS 

There are four types of boundary conditions to be considered. (1) 
specified temperature, (2) specified flux, (3) convection boundary layer, 
and (4) adiabatic. For a specified temperature boundary condition T^ is 
constant over boundary segment S^. I]or a specified flux q a "q over 
boundary segment S 2 . For a convection boundary layer q « h(T^ - T c ) over 
boundary segment Sg. Finally for an adiabatic boundary condition q = 0 
over boundary segment S 4 . 
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The boundary integral (2-21) becomes 



{£} "111 + Mfr} - {JTj 


( 2 - 22 ) 


hi - 


r \ 


(2-23) 


U 






¥\ 


kt [A]V/<<2 S>MdS 

vo 


(2-24) 


- 


kT 0 j [A]>>' ds 


(2-25) 



Si 



The integral on segment is zero because the variation of the functional 
was specified as zero on the boundary. The integral on segment S 4 is 
zero since there is no heat flow across the boundary. 

To obtain the extremum of the functional TT the nodal point tempera- 
tures ^T] which satisfy the following matrix equation must be found jjQ . 

(W-M)[t} +fc](t] - !%*} 0 (2-26) 



D. THE MATRIX EQUATION 

Let the time variable be ti such that t^ = iAt i = 0 , 1 , 2 ,.,. 
For this development i will refer to time increments. 



Let 



M - W-M 



(2-27) 



Equation (2-17) can be written as 

MW; .(0, 



(2-28) 
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where 



(2-29) 



Assuming that is constant over the interval £ t^t^ + i.e., 

itu - (m <41 

^ can be obtained by a Taylor^ expansion around t = t^ # 

■ {T};,, = It}; t At\T}. ♦- f\t}. (2-30) 

{T\ i* , « HI i ► Ai lT{ i + (2-31) 



frW- -\T] t +^(lT^, + fr}.) 



(2-32) 



Using equations (2-29) and (2-32) + ) and can be determined 

in terms of and • Solving (2-32) for gives 



itL, = r*(F3u. -UlO - It}, 



(2-33) 



Substituting this into (2-28) gives 



H+s&D l T k. = tflu) + CcO(Al'} : +fr] L ) 



(2-34) 



Also from equation (2-28) 



{tU = Ef’ffl - [cT'KKt}. 



(2-35) 



19 



Substituting (2-34) into (2-35) gives 






(2-36) 



or 




(2-37) 



In a multi-element body the element matrix equations must be assembled 
into a single matrix equation. The assembly and solution of the equation 
is a complex task but it can be greatly simplified by the use of a com- 
puter [6], 



E. FORMULATION OF THE ELEMENT MATRIX 

If the region to be analyzed is divided into plane triangular elements 
then a typical element would appear like the following figure. 




The temperature in an element can be described in terms of a linear com- 
bination of the coordinates of the elements nodal points and the values 
of the temperature at each of the nodal points. 

Assuming that the temperature distribution is of the form 

T(x,V,-L) - -Vo^x f4 3 C-Oy (2-38) 
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The temperature at each of the nodal points can be written as 

Tl - °l, + 

Tj - 0 (, + f V-o( 



J ^3 



(2-39) 



- <*, + x m <* t v Y m ^ 



Solving for the coefficients °( p 0 ^2 , anc * °^3 us *- n § Cramer's rule the 
following results are obtained. 



L err 



I Xi Vi 
I Xj Yj 
I "V 



- TUX 



(2-40) 



Note: 



is the area of the triangular element. 



<*. = 
°<z- 
°^3 “ 



T: 




yj 








Ivw ^ 


^ mr* 


/ 


it 




/ 




(i 


1 




VyV> 


1 


Kc 




) 






1 


X*. 


1 wv 



ZA 



l A 



ZA 
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I 



Vi 7; 

X *4 






Hvr\ 'Vvvi 

X,' v t - 



+“ ^W\ 




7 l 

7 - 













Recalling that the temperature distribution was previously described by 
the expression 



' ~ <mm (2-42) 

is the vector of the nodal temperatures and is the spacial 

dependence matrix which is a constant. Let 

(2-43) 

*.r °i/3 _ 

°(-LI ^3 

°l* °(J2. «3 2 

then 

T - <^>[a][t} 



<^> = <1 x y> 
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so 



^(% li+qfelj i 



T 



J 

\ 



^ + ^ L ) 
i^3, ^ ■+ vli + ^33 0 



> (2 



j 



T f <*iiTvJ + bu n +* lz jj +ctex,)x ^(^,»t+oi^7y ^iv J3 t^)v 

Therefore, comparing (2-43a) with (2-41) 



% - 


x j yj 
y„ 




V = 


x ^ ~y»? 

x; Y: 


= -li/* 


°l 1 3 ~ 


X; y £ 

*j i' 


= * *; - *,■ )i 




i z 

i Yj 


= ^ ^ 


°\Z7. 


i y i 

1 Y* 


- Y„ - Yj 


°*Z3 =~ 


1 $ 

1 V; 


- y, - y, 

u 'j 


^ 3( 


' *», 


r /U-*; 



-43a) 



23 



9f 



32. 



X; 



li - X 



Vw 



43 



/ 



XL 
A : 



Xj 



The spacial dependence matrix is then defined by 

(xyy^M;) ouy -xx^y/- yvo 



PQ = x 



2A 



JJ-Y. 


y - v- 


Y.-L' 






Xj-v : 



(2-44) 



PI 



F. THE METHOD OF LUMPING ELEMENT PROPERTIES 

When working with a computer it is usually more convenient to work 
with a quadrilateral element. The quadrilateral is divided into four 
triangles with the center nodal point common to each of the triangles, 
(see figure 2-2). 
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The matrix equation (2-28) can be assembled for this element in the 
following form 



K 11 K 12 K 13 K 14 K 15 




V 




Cn c 12 C 13 C 14 C 15 


k 21 K 22 K 23 K 24 K 25 




T 2 




C 21 C 22 C 23 C 24 C 25 


K 31 K 32 K 33 K 34 K 35 


4 


T 3 


> + 


C 31 C 32 C 33 C 34 C 35 


K 41 K 42 K 43 K 44 K 45 




T 4 




C 41 C 42 C 43 C 44 C 45 


^_5 1 K 52 K 53 K 54 K 55_ 




t 5 




C 51 c 52 c 53 c 54 c 55 ^ 


Partitioning this equation results in 



lT l 
To 



T. 






<; 3 >=s 






& 







> H- 


mm 


bj 


M 


^5j> K 55_ 


bj 




( C 5j> C 55_ 


I *5 


Uj 



i,j = 1,2, 3,4 

This equation can be multiplied out to yield two equations 



(2-45) 



[ K iIK T 4 + {***] X 5 + [ C iJ + M M (2 - 46) 

< K 5j>l T A + K 55 T 5 + < c 5j >^ii + C 55% - £ 5 (2-47) 

The specific heat matrix 0 =] is approximated by lumping the heat ca- 
pacities at the four external nodes. Thus DO becomes a diagonal matrix 
(4x4) and C 55 s± 0. Therefore (2-46) can be written as 



M W + { K i 5 } *5 + *{ £ i] ( 2 - 48 ) 

and (2-47) as 

<^5j) { T l] + K 55 T 5 - £ 5 (2-49) 

Solving (2-49) for T^ and substituting into (2-48) results in 

[ K J W + M k 55 ( f 5 - < K 5J>( T i} ) + [ C ij]l f i} -f £ i] < 2 ' 50) 
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or 



( Oy] - {\sl K 55 < K 5j) > { T i} + - f £ i] - { K i5} %5 £ 5 

Equation (2-51) is analogous to equation (2-28) except that jV] and 
are now (4x4) matrices and is a (4x4) vector |6^] . 



2-51) 

00 
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III. PROGRAM STRUCTURE 



This chapter describes the general composition of the program and 
gives a description of the function of each section of the program. 

A. PROGRAM IDENTIFICATION 

The program was written by Robert E. Nickell while in the Department of 
Civil Engineering, University of California at Berkeley. It was formerly 
identified as W2498,W035--2--57343, Nickell. It was written in July 1968 
and was revised by Allan B. Chaloupka in May 1969. Formulation was 
based upon a variational principle derived by M. E # Gurtin of Brown 
University B. 

The program has been renamed DFETHC-Finite Element Transient Heat 
Conduction. 

B. PURPOSE 

The purpose of the program is to produce high speed solutions to 
transient heat conduction problems. These problems may involve plane or 
axisymmetric geometry, non-homogeneous anisotropic material configura- 
tions, temperature dependent material properties, and time dependent 
boundary conditions. 

C. PROGRAMING INFORMATION 

The program was originally written in FORTRAN IV for the IBM 7094 
computer. It was converted for use on the IBM OS/360 Model 67 computer. 

The program was also rewritten in double precision (REAL *8). 

D. PROGRAM CAPACITY 

The size problem that can be analyzed by this program is subject to 
the following limitat ions’* 

Maximum number of nodal points 500 
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Maximum number of quadrilateral elements 



490 



Maximum number of materials ... .... 25 

Maximum difference of nodal point numbers 

for the same element «, . . . 31 



These are not stringent limitations and they may be adjusted some- 
what by the user to fit individual problems . The limits are changed by 
changing the appropriate DIMENSION and C0!°/ r N statements and one card in 
the segment that calculates the half -band width* The only real limita- 
tion placed upon the size of the problem that can be analyzed is the 
amount of core storage available in the computer used. 

E. GENERAL STRUCTURE 

The program consists of a main program and several subprograms used 
for repetitive operations. The main program can be divided into three 
parts. The first part is concerned with reading into the program and 
printing out the data describing each problem. It is portrayed by the 
flow diagram in figure 3-1. 

The second portion of the main program is concerned with the routing 
of operations according to certain logical variables. These variables 
describe four possible combinations. They are constant material proper- 
ties and constant boundary conditions, temperature dependent material 
properties and constant boundary conditions, constant material proper- 
ties and time dependent boundary conditions, and temperature dependent 
material properties and time dependent boundary conditions. Each of these 
possibilities is portrayed by a flow diagram in figures 3-2 through 3-5. 

The third portion of the main program is concerned with calculating 
the effective forcing vector, solving the matrix equation and printing 
the temperatures for each time increment. The operation is then directed 
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Figure 3-1 Main Program Flow Diagram, Part I 
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After the First 
Time Step r 



1 



Initialize Matrices 
and Vectors 




Figure 3-2 Main Program Flow Diagram, Part II 
Constant Properties and Boundary Conditions 
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Initialize Matrix 
and Vectors 



Form Heat Capacity 
and Conductivity Matrices 



Apply Constant 

Convection Boundary Conditions 



Apply Constant Temperature 
and Flux Boundary Conditions 



Form Effective 
Coefficient Matrix 



Figure 3-3 Main Program Flow Diagram, Part II 
Temperature Dependent Properties and 
Constant Boundary Conditions 



Put Matrix in 
Upper Triangular Form 
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Figure 3^4 Main Program Flow Diagram, Part II 
Constant Properties and Time Dependent Boundary Conditions 
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Figure 3-5 Main Program Flow Diagram, Part II 
Temperature Dependent Properties and 
Time Dependent Boundary Conditions 
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back to the beginning of the second section of the main program to start 
calculations for the next time increment. The third portion of the main 
program is portrayed by a flow diagram in figure 3-6. 

F. MAIN PROGRAM 

The main program is composed of a number of individual steps which 
comprise the finite element solution technique, and a system for the in- 
put and output of problem information. Not included in the main program 
are the matrix assembly processes, the matrix equation solving processes, 
and the boundary condition input routines,, These are included in sub- 
programs FORM, SYMSOL, TFUNC, CFUNC, FCBC, FTBC, and FFBC. 

1. Control Information 

The first information to be read into the program is the control 
information. This information sets the basic variables for each of the 
steps of the program and also routes the problem to various sections of 
the program. This information includes the number of nodal points, the 
number of elements, the number of convection boundary conditions, the num- 
ber of materials, the number of time sequences, whether the problem is in 
plane or axisymmetric geometry, the output print interval, the system of 
units, whether a final spacial temperature distribution is desired, and 
the reference temperature of the body. 

Also read in with the control information are two titles. The 
first is a 72 space title that can be used by the user to label each 
problem. The second is a 48 space title that is used to input the labels 
used with the system of units chosen. 

The appropriate titles and control information are printed on the 
first page of the printed output, for each problem. 
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Figure 3-6 Main Program Flow Diagram, Part III 
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2. Material Properties 



The thermal conductivity, density, specific heat, and heat gener- 
ation properties for each material are read into the program. The materi- 
al properties are then printed in order on the printed output page and 
labeled with the appropriate units. Since the program is for anisotropic 
materials, k^jkyy, and are read in order to form the thermal con- 

ductivity tensor. The specific heat and density are read in as a product 
in order to reduce the number of variables in the program. 

3. Time Sequencing Information 

In this section the time sequencing information is read in for 
each sequencing scheme. This information consists of two logical varia- 
bles and two numerical variables. The logical variables indicate whether 
there are temperature dependent properties and if there are time dependent 
boundary conditions. The two remaining variables describe the number 
of time steps and the size of each time increment. The time sequencing 
information is printed on the printed output page with appropriate labels. 

4. Nodal Point Information 

Ordinarily this segment of the program will read the x,y coordinates 
and initial temperature of each nodal point and print the information 
with the appropriate dimensions. Because of the possibility of a large 
number of nodal points (500) and the equally large possibility of human 
error in punching data cards there is a segment of the program that will 
generate nodal point information. This segment has a few limitations, 
which will be discussed in chapter IV, but it is usually able to generate 
almost all of the interior nodal points. As an example in the configura- 
tion shown in figure 3-7 the nodal points indicated by squares were read 
into the program and the remaining nodal points were generated by the 
program. 
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Figure 3-7 Nodal Point and Element Generation 
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5. Element Information 



The purpose of this segment is about the same as the previous 
one. The number of the element, the number of each of its nodes, the 
type of material, and the amount of heat generation within the element 
are read into the program and printed out with appropriate labels. Again 
this segment has the capability to facilitate the loading of data by 
generating element information. This segment also has some limitations 
which will be discussed in chapter IV, but it is usually able to generate 
most of the elements in a given configuration. As an example, in figure 
3-7 all the elements were generated except those that are starred. 

Also included in this segment is a routine that calculates the 
half -band width of the problem. This is an important control parameter 
that is used later in the program to assemble the matrix equation. 

6. Time Steps 

The time step portion of the program includes two large loops. 

The first loop executes the program through each time sequence of the 
problem. The second loop executes the program through the required num- 
ber of time increments for each time sequence. Inside these two loops 
the two logical time sequencing variables and the loop counters control 
the routing of the problem. 

The two logical variables control the choice of variable material 
properties and boundary conditions. If the material properties and bound- 
ary conditions are constant then the coefficient matrix and forcing vector 
are calculated only once. The matrix equation is then solved for suc- 
cessive time increments using only the original matrices. If variable 
material properties are involved the effective coefficient matrix must 
be reapplied since the coefficient matrix is reformed from zero values. 

If time dependent boundary conditions are involved then the contribution 
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of the boundary conditions to the coefficient matrix must be reapplied. 
The coefficient matrix from the original material properties is kept in 
temporary storage. When both time dependent boundary conditions and 
temperature dependent properties are used then the coefficient matrix 
and forcing vectors must be reformed for each time increment based upon 
the immediate time and temperature values and initial properties. The 
data for temperature dependent properties is calculated from the initial 
conditions and the present nodal temperatures. The data for time depend- 
ent boundary conditions is either read into the program with each time 
increment or generated within the program. 

Different time sequences can be run for the same problem. This 
enables the user to speed or slow the time change during a particular 
segment of the problem. 

7. Convection Boundary Conditions 

A convection boundary condition occurs when the problem includes 
a convection boundary layer as part of its boundary conditions. One con- 
vection boundary condition consists of two nodal point numbers, the con- 
vective heat transfer coefficient and the ambient fluid temperature. The 
two nodal points define the convection boundary for each condition. 

Time varying convection boundary conditions can be obtained by 
varying the ambient fluid temperature. The new fluid temperature may be 
read with each new time step or it may be generated by using the function 
subprogram FCBC. FCBC is used in the same manner as the function sub- 
programs for temperature dependent material properties. The argument of 
FCBC is TIME and the value of FCBC is currently set at one. If tempera- 
ture dependent material properties are used then the convection boundary 
condition section stores and reapplies the convection boundary condition 
contributions to the coefficient matrix and the forcing vector. 
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8 . Temperature and Flux Boundary Conditions 

A temperature or flux boundary condition occurs when the problem 
specifies a temperature or flux at a nodal point. One temperature or flux 
boundary conditions consists of the nodal point number and the value of 
the temperature or flux at that nodal point. 

Time varying temperature and flux boundary conditions can be ob- 
tained by varying the temperature or flux at each nodal point. The new 
temperature or flux may be read with each new time step or it may be 
generated by using the function subprograms FTBC or FFBC for time depend- 
ent temperature or flux respectively. FTBC and FFBC are used in the 
same manner as the function subprograms for temperature dependent materi- 
al properties. The argument of each is TIME. Their current values are 
both one. If temperature dependent properties are used the temperature 
and flux boundary condition section stores and reapplies the boundary 
conditions to the coefficient matrix and the forcing vector. 

9. Temperature Output 

This segment of the program is concerned with the print out of 
the temperature solutions for each nodal point at each time step. The 
value of the time is printed and the nodal point numbers and correspond- 
ing temperatures for that time are printed in rows below. If it is de- 
sired, the temperature can be printed with spacial data at the end of 
each time sequence. 

G. SUBROUTINE FORM 

The subroutine FORM is used to calculate the effective element ther- 
mal conductivity, heat capacity, and heat generation. The subroutine 
performs these operations for each of the elements of the configuration 
and, if temperature dependent properties are used, for each time increment. 
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Figure 3-8 The Matrix Equation 




Figure 3-9 Computer Storage for the 
Matrix of Coefficients 
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First the mean x,y coordinates and the mean temperature of the element 
are calculated by averaging the values of each of the nodes. Using the 
mean x,y coordinates the quadrilateral element is divided into four tri- 
angles. The conductivity tensor, heat capacity, and heat generation 
vectors are formed for each of the four triangular sections. As these 
matrices are formed they are assembled into the larger (5) nodal point 
equation for the quadrilateral element. The heat capacity and heat gen- 
eration vectors are assembled so that the fifth term of each vector is 
dropped when forming the lumped properties matrix equation. The center 
nodal point is eliminated from the thermal conductivity tensor (5x5), 
according to equation (2-52) , to give the thermal conductivity tensor 
(4x4) for the lumped properties matrix equation. The thermal conductivity 
tensor and heat capacity and heat generation vectors for each element 
are returned to the main program where they are loaded into the coeffi- 
cient matrix and vectors used in the solution of the total matrix equa- 
tion for each time step. For problems where the properties of the 
materials are not temperature dependent the thermal conductivity tensor 
and the heat capacity and heat generation vectors are calculated only at 
the beginning of the problem. If temperature dependent material proper- 
ties are included these matrices are calculated at the beginning of each 
time step. 

H. SUBROUTINE SYMSOL 

The subroutine SYMSOL in conjunction with the loading methods of the 
main program is used to solve the finite element matrix equation 0 It 
consists of two parts. The first part reduces the coefficient matrix 
to its final form. The second part forms the effective forcing vector 
and solves for the temperatures at each of the nodal points. 
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In an ordinary formulation the final matrix equation for the finite 
element method would look like figure 3-8. The matrix of coefficients 



core storage, the matrix of coefficients is stored in columns. The first 
column is the diagonal of the matrix of coefficients. (See figure 3-9). 
All the operations performed on the matrix of coefficients are done on 
the upper half of the symmetric banded matrix in its stored form. 

After the coefficient matrix is formed it is put into upper triangular 
form. This yields a matrix equation of the form shown below. 



This is accomplished by the first segment of SYMSOL, The second portion 
of SYMSOL, using the effective forcing vector, starts with the temperature 
at the bottom of the nodal temperature vector and back substitutes to 
solve for the nodal temperatures. These operations are performed for each 
time increment. 




banded symmetric matrix. In the computer, in order to conserve 




Figure 3-10 Matrix Equation in Upper Triangular Form 
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I. TEMPERATURE DEPENDENT PROPERTIES 



There are two function subprograms that are used in conjunction with 
temperature dependent properties. They are TFUNC and CFUNC. TFUNC is 
used for the temperature dependent thermal conductivity and CFUNC is used 
for the temperature dependent product of the density and specific heat. 

In the programs present configuration each of these function subprograms 
is prepared for temperature independent properties and their values are 
therefore equal to one. If the thermal conductivity were temperature 
dependent in the form of: k = k Q (l + bT) , where b is a constant coeffi- 
cient and k Q is the value of the thermal conductivity at T * 0, then 
TFUNC would have to be modified by the user to include the factor 1 + bT. 
The material properties read into the program would then be zero tempera- 
ture properties. 
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IV. USER INFORMATION 



This chapter describes the manner in which problems should be formu- 
lated. Also the procedures used in inputting data into the program are 
inc luded. 

A. INPUT DATA 

Data for each problem processed by this program is assembled in groups 
The first group includes the titles and all the control parameters,, The 
second group is the material properties. The third group is the time se- 
quencing information. The fourth group is nodal point information. The 
fifth group is element information. The sixth group is boundary condi- 
tions. After each data deck there must be three blank cards. 

B. CONTROL PARAMETERS 

The following control parameters are read into the program. 



NUMNP 


number of nodal point 


NUMEL 


number of elements 


NUMCBC 


number of convection boundary conditions 


NUMMAT 


number of materials 


NSEQ 


number of time sequences 


KAT 


type of geometry, 0 implies plane geometry 


INTER 


output print interval 


UNIT 


units system 


KPUNCH 


spacial temperature output 


JUMP 


generated temperature of flux boundary conditions 


KEY 


generated convection boundary conditions 


TO 


initial temperature 
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Also included with the control data are two labels, TITLE (18) is an 



arbitrary label that may be used to identify individual problems. It is 
72 spaces long. The second label is UN(12). UN(12) is an array of four 
unit labels. Each has 12 spaces. The units must be printed in the fol- 
lowing order: length, mass, time, and temperature. These labels must be 
left justified. 

C. MATERIAL PROPERTIES 

The following material properties are read into the program. 



XCOND 


k 




XX 


YCOND 


kyy 


XYCOND 


k 

xy 


SPHT 


e c 


HX 


heat generation 



Material properties are read in with one material to a card. If the 
material is isotropic k xx = kyy =0. If the material is anisotropic the 
three values will not be the same. Each material is given a number ac- 
cording to the order it is read into the program. 

If temperature dependent material properties are used both logical 
variables must be true. Also both JUMP and KEY must be some number other 
than zero and FCBC,FTBC, and FFBC must be one. 

D. TIME SEQUENCING INFORMATION 

The following time sequencing information is read into the program. 



LTAG 


temperature dependent properties (true or false) 


MTAG 


time dependent boundary conditions (true or false) 


ITAG 


number of time steps 


TAG 


time increment 



The time sequencing information controls the operation of the program 
through the use of the logical variables LTAG and MTAG but it can also 
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be used to slow or speed time during the problem. An example would be 
two sequences with the first solving for the temperature 60 times during 
the first hour of the problem (i.e., every minute). The second time se- 
quence would solve for the temperature 6 times in the second hour of prob- 
lem time (i.e., every 10 minutes), 

E. NODAL POINT CONSTRUCTION 

The following nodal point information is read into the program for 
nodal point data. 

N number of the nodal point 

KODE 1 indicates a temperature boundary condition 

X x coordinate of the nodal point 

Y y coordinate of the nodal point 

T initial temperature of the nodal point 

The data for each nodal point is put on one card. 

Earlier it was mentioned that the program could generate nodal points. 
This feature has a few limitations that must be observed or incorrect 
information will result. A nodal point whose KODE is 1 cannot be gener- 
ated. Nodal points must be generated along straight lines. The initial 
temperature of the generated nodal points must be the same and equal to TO, 
the initial temperature. To generate nodal points simply leave out the 
data cards between two nodal points on a straight line and the program 
will generate the information for the nodal points in between. 

In constructing the mesh of nodal points for a problem it is sometimes 
impossible to divide the area to be analyzed into uniform rectangular 
quadrilaterals. At times it will not be desirable to do so. Consider 
figure 4-1. In preparing a nodal mesh for this configuration the first 
step would be to assume that by symmetry there is no heat flux across the 
center lines. This reduces the figure to one of its quadrants as in 
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Figure 4-1 Irregularly Bounded Configuration 
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Figure 4-2 Reduced Irregularly Bounded Configuration 



figure 4-2. The size of the quadrilateral elements would vary depending 
upon how interested the user was in a certain area of the figure. In 
figure 4-2 the area around the curve is being more closely examined but 
the mesh size must also be small to enable the user to approximate any 
curved surface using a series of straight line segments. It is almost 
impossible to divide an area with a curved boundary without encountering 
a right triangular element. To convert this to a quadrilateral element 
place a nodal point in the middle of the hypotenuse. 




Figure 4-3 Conversion of a Right -Triangle 
to a Quadrilateral 

In numbering the nodal points it must be kept in mind that the maxi- 
mum difference in nodal point numbers in one element is 31. This is a 
limitation set by the ha If -band width of the computer program. It is 
convenient to increase the numbering of the nodal points and the elements 
in the same pattern, although it is not a requirement. 

F. ELEMENT CONSTRUCTION 

The following element information is read into the program for ele- 
ment data. 

N number of the element 

IX(N,K) K *= 1 - 4 numbers of the element f s nodal points 
IX(N,5) material type 

HTGEN element heat generation 

The data for each element must be placed on one card. The nodal point 
numbers are placed in clockwise order starting with the number in the 
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upper left hand corner. Elements may be generated by the program but 
there are some limitations. Only elements of the same material type can 
be generated at one time. The heat generation must be the same for the 
elements generated together. Elements must be generated in rows or col- 
umns, To generate elements simply leave out those element cards that 
you wish generated. 

G. BOUNDARY CONDITIONS 

There are two types of boundary conditions that may be read into the 
program, convection and temperature or flux. (An adiabatic boundary con- 
dition is a flux boundary condition with the flux equal to zero.) The 
following information is read into the program for convection boundary 
conditions. 

IB nodal point number 

JB nodal point number 

HB convective heat transfer coefficient 

TEMPB ambient fluid temperature 

The data for each portion of the convection boundary is put one one card. 
If time varying ambient temperature is to be generated by the program the 
control parameter KEY must be some number other than zero, say one. If 
the time varying ambient temperature is not generated by the program new 
data for each segment of the convective boundary must be read into the 
program with each new time step. In this case KEY is zero. 

The following information is read into the program for temperature 
and flux boundary conditions. 

NN nodal point number 

TQ temperature or flux 

The data for the temperature or flux boundary condition at each nodal 
point is put on one card. If time varying temperature or flux is to be 



52 



generated by the program the control parameter * JUMP must be some number 
other than zero. If time varying temperature or flux is not generated! 
by the program new data for each boundary nodal point must be read into 
the program with each new time step* In this case JUMP is zero. Both 
temperature and flux boundary conditions must be introduced into the pro- 
gram in the same manner. For both convection and temperature or. flux 
boundary conditions, if time variations are used the logical variable 
MT AG must be true. Also, regardless of the position of the last nodal 
point, a boundary condition card must be read into the program for it. 

If the last nodal point is not a boundary point then the data card con- 
sists of the nodal point number and a value of zero in the temperature 
or flux field. 

H # . UNITS 

The units used for each problem must be consistant with the input 
data or the labels on the output will be incorrect. Presently there are 
six different units systems available. Each system consists of a length, 
mass, time, and temperature. Each system has a number which must be read 
in as a control parameter (UNIT), The systems are: 



1. 


feet, pound mass. 


hours, °F 


2. 


inches, pound mass 


, seconds, °F 


3. 


feet, pound mass. 


seconds, °F 


4. 


meters, kilograms. 


seconds, °C 


5. 


centimeters , grams 


, seconds, °C 


6. 


feet, pound mass, 


minutes, °F 



If the user does not prefer the units available he can use another 
system in his input and change the labels or just ignore the labels al- 
together. The system of units will depend greatly upon the configuration 
of the problem and the temperature differences involved. 
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I. ERROR EXITS 



There are three error exits in the program. The first occurs if a 
non-positive number of nodal points is read into the program. This exit 
is used to stop the program. The other two occur in the nodal point and 
element information reading sections. If these statements are executed 
then the execution of the program is terminated. If either of the last 
two exits are executed a statement will be printed to indicate which one 
it was. Usually termination means a bad data card or that the user was 
trying to generate points beyond the capability of the program. 

J. COORDINATES 

The cartesian system of coordinates x and y are used for problems in- 
volving plane geometry. When axisymmetric problems are considered the 
x, O ', 2 system of polar coordinates is used. The analysis of axisymmetric 
problems is conducted in the x,z plane. 

K. TEMPERATURE OUTPUT 

The temperature output can be modified in two ways. The number of 
time increments printed out can be limited by using the control parameter 
INTER. If the program is directed to solve for the temperature every 
minute for 100 minutes and INTER is equal to 25, temperature information 
will only be printed for 25, 50, 75, and 100 minutes. The second manner 
in which the output may be modified is at the end of each time sequence 
the x,y coordinates of each nodal point and the temperature are printed 
if the control parameter KPUNCH equals one. 

L. INPUT DATA PREPARATION 

The following sections describe the manner in which the input data 
is to be placed on data cards. The data is prepared in groups and the 
final data deck is assembled by putting the data groups together in order 
of their listing below. 
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1. Control Information 



(a) Title Card (18M) : 72 arbitrary spaces for an alphameric 

problem label punched in columns 1-72 

(b) Units Card (12M): Alphameric unit labels 



Columns 


1-12 


length 






13-24 


mass 






25-36 


time 






37-48 


temperature 


Labels should 


be left justified. 




(c) Control 


Information 


Card (1115, F10.0) 


Co lumns 


1-5 


NUMNP 


number of nodal points 




6-10 


NTJMEL 


number of elements - 




11-15 


NUMCBC 


number of convection 








boundary conditions 




16-20 


NUMMAI 


number of materials 




21-25 


NSEQ 


number of time sequence 




26-30 


KAT 


type of geometry 




31-35 


INTER 


print interval 




36-40 


UNIT 


units system (integer) 




41-45 


KPUNCH 


spacial distribution 




46-50 


JUMP 


temperature conditions 




51-55 


KEY 


convection conditions 




56-65 


TO 


initial temperature 



"I" formats must be right justified. 

2. Material Properties 

(a) Material Properties Card (5D10.3) 
Columns 1-10 XCOND kxx 
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11-20 



YCOND 



k yy 

21-30 XYCOND k xy 

31-40 SPHT pc 

41-50 HX heat generation 



3 . 



One card for each material. 

Time Sequencing Information 

(a) Time Sequencing Card (2L5, I10,D10. 3) : 



Columns 1-5 


LT AG 


temperature dependent 






properties (true or false) 


6-10 


MTAG 


time dependent boundary 






conditions (true or false) 


11-20 


ITAG 


number of time steps 


21-21 


TAG 


time increment 



4. 



5 . 



One card for each time sequence. 
Nodal Point Information 
(a) Nodal Point Card (2I5,3F10.0) 
Columns 1-5 N 

6-10 KODE 

11-20 X 
21-30 Y 
31-40 T 

One card for each nodal point not 
E lement Information 
(a) Element Card (615,010.3): 
Columns 1-5 N 

6-10 IX(N,1) 



nodal point number 
temperature boundary 
condition 
x coordinate 
y coordinate 
initial temperature 
generated by the program. 



element number 
nodal point number 



56 



11-15 



IX (N, 2) nodal point number 

16-20 IX(N,3) nodal point number 

21-25 IX (N, 4) nodal point number 

26-30 IX(N,5) material number 

31-40 HTGEN element heat generation 

One card for each element not generated. Nodal points are num- 
bered clockwise starting in the upper left hand corner. 

6 . Convection Boundary Conditions 

(a) Convection Boundary Card (2I5,D10. 3,F10. 0) 

Columns 1-5 IB nodal point number 

6-10 JB 

HB 



11-20 



21-30 



nodal point number 
convective heat transfer 
coefficient 

ambient fluid temperature 



TEMPB 

One card for each convection boundary condition. 

7. Temperature and Flux Boundary Conditions 

(a) Temperature or Flux Boundary Card (I10,D10.3) 

Columns 1-10 NN nodal point number 

11-20 TQ temperature or flux 

One card for each temperature or flux boundary point. A tempera- 
ture or flux boundary card for the last nodal point must be read even 
if it isn f t a boundary point. 

8. Multiple Problems 

The above cards comprise a data deck for one problem. Multiple 
problems may be solved by putting their data decks together. The only 
limitation on this procedure is computer time. It is important that 
there be three blank cards after the last problem data set in the data 
deck. This is to terminate operation of the program. 
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RECOMMENDATIONS 



Since the number of problem types that can be solved by this program 
is limited mostly by the imagination of the user it is highly recommended 
that this engineering tool should be brought to the attention of those 
students studying in the area of unsteady heat conduction. Only through 
complete familiarization with the program can the user make use of its 
great potential to solve problems involving complex configurations and 
boundary conditions. 

Recommendations for further study in this area are: 

1. Expansion of the finite element method to three dimensional con- 
figurations. 

2. Revision of the boundary condition segments of the program to 
make use of random access disks for temporary storage. 

3. Construction of a supplementary program that will punch input 
data cards for complex time dependent boundary conditions. 

4. An isothermal contour plotting subroutine capable of being used 
in conjunction with the present NPS Computer Center subroutine DRAW. 

5. Use of the program in conjunction with the NPS Computer Center 
IBM Optical Display Unit, in order to observe transient heat conduction 
differences after minor alterations in mesh configuration. 

6. Use of the program in conjunction with existing finite element 
programs in the area of stress analysis to obtain thermal stresses. 
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APPENDIX A 



In this section some solutions to different configurations will be 
presented as examples of possible problem types that may be solved with 
this program. Also the results of these examples are compared with 
analytic solutions of the same configurations. In order to avoid con- 
fusion the same material nickel (99.9%) and the same initial temperature, 
100°F, are used in each problem. 

The material properties of nickel are: 

^ = 556 lbm/ft 8 
c = . 1065 BTU/lbm-°F 
- 882 ft 2 /hr 

Example 1. 

This is a problem in one dimensional heat conduction. It could also 
be considered a problem of heat conduction in a semi-infinite plate by 
assuming that lines parallel to the direction of heat flux are adiabatic 
boundaries. Because of this assumption we only need one row of elements 
across the plate which is one foot wide (see figure A-l), The top and 
bottom sides of the elements are adiabatic boundaries. Also the end 
boundary conditions are adiabatic on the left and are constant 0°F on the 
right. The nodal points are numbered from left to right, starting with 

the top row. The elements are also numbered from left to right. The 
computer solution provides answers at distinct time intervals while most 




Figure A-l One Dimensional Heat Conduction 
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analytic solutions are given in terms of the dimensionless time parameter 

2 

Fourier number, NpQ. The Fourier number, N^q, equals 0^ t/L , where T 
is the thermal diffusivity of the material. The comparisons will be made 



at various Fourier numbers. The analytic solution was obtained from 
page 98 of reference &]. 



Nodal 


NFO 


Cars law and 


Finite 


Point 




Jaeger 


Element 


1 


.04 


100 


99.6 


5 


.04 


99 


99.0 


9 


.04 


96 


95.9 


13 


.04 


85 


84.7 


17 


.04 


54 


54.5 


1 


. 1 


95 


94.8 


5 


.1 


93 


91.5 


9 


.1 


82 


82.5 


13 


. 1 


63 


63.6 


17 


. 1 


35 


35. 1 


1 


.4 


48 


48. 1 


5 


.4 


45 


45.8 


9 


.4 


38 


38.9 


13 


.4 


28 


28.3 


17 


.4 


15 


14.9 



Example 2. 

This example has the same mesh configuration and boundary conditions 



as the first example. The initial temperature distribution is changed 



to the 


form T a 100 - 50x°F. 


The analytic solution was 


obtained from 


page 98 


of reference QT). 






Nodal 


%0 


Carslaw and 


Finite 


Point 


Jaeger 


Element 


1 


.01 


95 


94.0 


5 


.01 


90 


90.0 


9 


.01 


80 


79.8 


13 


.01 


70 


68.0 


17 


.01 


60 


53.0 


1 


.04 


99 


89.0 


5 


.04 


86 


85.6 


9 


.04 


77 


76.9 



60 



13 


.04 


62 


62 . 1 


17 


.04 


36 


37.0 


1 


. 1 


80 


80.0 


5 


. 1 


76 


76.0 


9 


. 1 


67 


66.3 


13 


. 1 


50 


49.7 


17 


. 1 


28 


27.0 



Example 3. 

This example is a solution for two-dimensional heat conduction. The 
configuration is a square plate with an initial temperature of 100°F and 
boundary conditions of 0°F on all sides. Using symmetry the user need 
only examine one fourth of the plate, i.e. one corner (see figure A- 2). 




Figure A-2 Two Dimensional Heat Conduction 



The nodal points are numbered from left to right and from the bottom row 
to the top. The elements are numbered in the same manner. In this case 
the problem consists of a plate with two sides at 0°F and the other two 
insulated for an adiabatic boundary condition. The analytic solution is 
obtained by using the solution for one dimensional heat conduction in a 
semi-infinite plate and the method of multiplicative superposition. 



61 



The analytic solution was obtained from page 118 of reference Qf] . 



Nodal 


Npo 




Finite 


Point 




Chapman 


Element 


8 


.01 


27.0 


30.0 


15 


.01 


70.6 


70.6 


22 


.01 


92.0 


90.5 


29 


.01 


97.0 


97.3 


36 


.01 


100.0 


98.8 


8 


.05 


5.75 


6.29 


15 


.05 


21.2 


22.2 


22 


.05 


40.4 


40.8 


29 


.05 


54.0 


55. 1 


36 


.05 


59.3 


60.4 


8 


.1 


1.96 


2.23 


15 


.1 


7.3 


8.1 


22 


.1 


14.4 


15.2 


29 


.1 


19.4 


21.1 


39 


.1 


22.1 


23.3 


Example 4. 









In this problem the same mesh configuration as the previous example 



is used. The boundary conditions are changed so that the sides are adia- 
batic, the top is at a constant 100°F. The initial temperature is 100°F. 
By taking an extremely large time increment, that is of n magnitude 
greater than the maximum number of significant figures available in the 
computer, say 10^ f the first time step will yield t. ho steady state 
solution to the problem. This is made possible because the I ime depend- 
ent contributions to the coefficient matrix and the forcing vector are 
multiplied by the inverse of the time increment which in this case is 
essentially zero (10 ^ ). The solution to this problem Is a linear 
temperature distribution and the results below were obtained by using a 
time Increment of 10^ minutes. 
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Point Element. 
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APPENDIX B 



A 


coefficient matrix 


ACALTR 


dummy logical variable 


B 


effective forcing vector 


CONDI 


k x in an element 


CONDJ 


ky in an element 


CONDJ 


k X y in an element 


H 


convection boundary layer heat transfer coefficient 


HX 


heat generation in a material 


INTER 


output print interval 


ITAG 


number of time steps 


IX 


k ^ 1 - 4 nodal point numbers 



k s 5 material number 



JUMP 


temperature and flux boundary generation 


KAT 


geometry type 


KEY 


convection boundary generation 


KODE 


temperature or flux boundary parameter 


KPUNCH 


spacial temperature distribution parameter 


LPRINT 


print parameter 


LTAG 


logical variable 


MB AND 


half band width 


MTAG 


logical variable 


MTYPE 


material number 


NSEQ 


number of time sequences 


NUMCBC 


number of convection boundary conditions 


NUMEL 


number of elements 


NUMMAT 


number of materials 
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NUMNP 

SPHT 

T 

TAG 

TEMP 

TIME 

TITLE 

TMEAN 

TO 

UN 

UNIT 

X 

XCOND 

XMEAN 

XYCOND 

Y 

YCOND 

YMEAN 



number of nodal points 

the product 

nodal temperature 

time increment 

ambient fluid temperature 

time variable 

identification title 

mean temperature in an element 

initial temperature 

unit labels 

unit system parameter 

nodal point x coordinate 

kjj in a material 

mean x coordinate of an element 

kxy in a material 

nodal point y coordinate 

ky in a material 

mean y coordinate of an element 
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APPENDIX C 



Title and Labels 

THESIS CHECK PROBLEM HEAT CONDUCTION IN A RECTANGLE 
FEET POUND MASS MINUTES FAHRENHEIT 



* 

Control Information 

36 25 C 1 1 C 1 6 1 0 0100.0 



* 

Material Properties 

+8.660D-01+8.660D-01+0. 000D+00+5. 930D+0 1+3.000D+00 
* 

Time Sequencing Information 
FALSEFALSF 250+ 1 • OCCC-C 1 



Nodal 


Point 


Information 




1 


1C.0 


0.0 


100.0 


2 


10.1 


0.0 


100.0 


•a 


10.2 


C.O 


100.0 


4 


10.3 


0.0 


100.0 


5 


10.4 


0.0 


100.0 


6 


10.5 


C.O 


100.0 


7 


10.0 


0 . 1 


100.0 


12 


00.5 


0.1 


100.0 


13 


10.0 


0.2 


100.0 


18 


00.5 


0.2 


100.0 


19 


10.0 


0.3 


100.0 


24 


00.5 


0.3 


100.0 


25 


10.0 


0.4 


100.0 


30 


00.5 


0.4 


100.0 


31 


10.0 


C.5 


100.0 


36 


00.5 


0.5 


100.0 



Element Information 



1 


7 


8 


2 


1 


1+0. OOOD+OO 


5 


11 


12 


6 


5 


1+0. 0000+00 


A 


13 


1^ 


8 


7 


1+0.0000+0C 


10 


17 


18 


12 


11 


1+0. 0000+00 


11 


19 


20 


14 


13 


1+0. 0000+00 


15 


23 


24 


18 


17 


1+0. 0000+00 


16 


25 


26 


20 


19 


1+0.0000+00 


20 


29 


30 


24 


23 


1+0.0000+00 


21 


31 


32 


26 


25 


1+C. 0000+00 


25 


35 


36 


30 


29 


1 + 0. 0000+00 



% 

Boundary Conditions 

1 + 0. 0CCD + 00 
2+C. 0CCD+00 
3+C. 00CD+00 
4+0. 00CD+00 
5+0.0000+00 
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6+C.C0CC+00 

7+c. occr+oo 

13+C.OCGD+CO 

i°+o.occn+oo 

25 + C. OCCD+OO 
31 + 0. CCCP+CC 
36+0. 0CCC+00 



NOTE: Starred Lines are not Data Cards. 
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